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A FOURTH-ORDER BOX METHOD FOR SOLVING 
THE BOUNDARY-LAYER EQUATIONS 

By Stephen F. Wornom 

Langley Research Center 


SUMMARY 

A fourth-order box method for calculating high-accuracy numerical 
solutions to parabolic, partial differential equations in two variables 
or ordinary differential equations is presented. The method is the natural 
extension of the second-order Keller "Box" scheme to fourth order and is 
demonstrated with application to the incompressible, laminar and turbulent 
boundary-layer equations. Numerical results for high-accuracy test cases 
show the method to be significantly faster than other higher-order and 
second-order methods available today. 

INTRODUCTION 

Recently, attention has been given to numerical methods for solving 
the boundary-layer equations which have truncation errors that are of 
higher order than the second-order methods presently in use today. The 
term higher order will refer in this paper to the truncation error in the 
coordinate normal to the body surface. The truncation error in the tan- 
gential coordinate is second order. The advantages of higher-order methods 
are twofold. P irst, they can be used to obtain a numerical solution as 
accurate as a second-order method with considerably less computer time and 
storage or, alternatively, they may be used to produce a significantly more 
accurate solution for the same amount of run time and storage as a second- 
order method. 

The higher-order methods proposed thus far for the boundary-layer equa- 
tions consist of three-point schemes which fall into two classes. The first 


class consists of collocation procedures which are fourth order for uniform 
meshes. These procedures treat the functional value and certain derivatives 
as unknowns at three collocation points and can be derived via Taylor Series 
(Hermite) or polynomial interpolation (Spline). In this category are the 

I 1 O O 

Pade approximation of Kreiss or so-called compact scheme , the Mehrstellen 0 
formulation, and the formulation of Peters^. The second class consists of 
methods for variable meshes. In this category are the Spline collocation 
methods of Rubin and Graves 5 and Rubin and Khosla 5 *^. Both classes are 
similar in that the resulting finite-difference equations involve three 
nodal points, but are different in that the first class is restricted to 
constant meshes whereas the second class is applicable to variable meshes. 

One disadvantage of the higher-order methods involving three nodal 
points is that the usual boundary conditions for incompressible flow, 
u = v = 0 at the surface and u + u g as the boundary layer merges with the 
the mainstream, are no longer sufficient since the resulting system of 
finite-difference equations contains two more unknowns than equations. 

This difficulty is usually finessed by adding an additional boundary 
condition at the cuter edge of the layer and an additional equation or 
boundary condition at the surface boundary. The choice of which additional 
boundary conditions or equations to use is not clear. 

Another disadvantage that is evident in some higher-order methods is 
the complexity of the resulting nonlinear finite-difference equations and 
the associated difficulty of solving them efficiently at each column. For 
example, the "Spline 4" method of Rubin and Khosla 5 ^ would seem to require 
the solution of a 5X5-block tridiagonal matrix in order to solve the fully- 
coupled incompressible boundary-layer equations. A simpler solution scheme, 
which "lags" the continuity equation one iteration behind the momentum equa- 
tion, is reported by Rubin and Khosla 6 . Since the equations in this scheme 
are solved uncoupled, the errors will diminish in a linear manner with each 

iteration. In contrast, the better second-order methods, such as the Davis 

ft 9 * 

Coupled Scheme (DCS) or the Keller Box Scheme (KBS), solve the equations 

fully coupled with Newton iteration and thus, for laminar flows, quadratic 

convergence is achieved. Hence, the advantages of some higher-order methods 

relative to second-order methods, may be diminished or lost entirely in 

practical engineering calculations. 


The purpose of this report is to present a fourth-order box method for 
obtaining numerical solutions to parabolic, partial -differential equations in 
two variables or ordinary differential equations, with application here to 
the steady, two-dimensional, incompressible, laminar or turbulent boundary- 
layer equations. The method has the following features: (1) The method 

results in finite-difference equations that involve only two nodal points 
and therefore is formally fourth-order accurate on all grids, (2) the method 
results in a 3X3 matrix of unknowns at each nodal point when the equations are 
solved in a coupled manner, (3) the method utilizes Newton iteration and 
demonstrates quadratic convergence for laminar flows, and (4) the method 
requires only the standard boundary conditions u = v = 0 at the body sur- 
face and u + u. as the boundary layer merges with the mainstream. In short, 

c 

the method is the natural fourth-order extension of the second-order Keller 
Box Scheme; it is an example of a general class of high-accuracy, two-point 
methods discussed briefly by Keller*^. Keller offers an operation-count 
analysis that suggests that such methods may be superior to the well-known 
KBS (with Richardson extrapolation) to achieve accuracy less than or equal 
to 0(h®). The validity of this conjecture is borne out in the present 
formulation and numerical results. 
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SYMBOLS 

3X3 matrices defined by equation (27) 

i ,jth element of the A n 

i,jth element of the B n 

3X3 matrix given by equation (31) 
surface skin friction coefficient, 

c f - C/(v2 °-V 2 ) 

value of C f for 640 intervals across the boundary 
layer 
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constants in grid stretch function - see equation (42) 
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s exact 


3X3 matrix given by equation (37) 

coefficients in the solution technique - see 
equations (32) to (34) 

percentage error in the wall shear - see equation (41) 

damping factor - see equation (10c) 

normalized longitudinal velocity component in the 
boundary layer f * u/U e 

represents any sufficiently differentiable quantity 
step size in the n coordinate h = An ^ = n n - n n _-| 

von Karman constant, k = .41 
reference length 
l = 1 + £ 

components of the vector - see equation (26) 
characteristic Reynolds number R = p u L /u 

6 oo ct' oo 

af/an 

value of s-| (wall shear) for 640 intervals across the 
boundary layer 
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defined by equation (10b) 
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u 
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Y + 
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z n 
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nondimensional inviscid flow velocity in the x-direction 
* . * 

nondimensional time-averaged viscous flow component in 

ft ft 

the x-direction u - if /U m 
friction velocity u^’ ' = VlV I A- 

nondimensional time-averaged viscous flow component in 
the Y-direction 7 v*/U w * 

transformed viscous flow component in the y-direction - 
see equation (7) 

vector defined by equation (26) 

nondimensional position coordinate measured along body 

ft ft 

surface from leading edge or stagnation point, x - x /L 
stretched normal coordinate Y y 

law of the wall coordinate Y + - u t *y*/v co 

nondimensional physical distance normal to body surface 
* * 
y = y /l 

vector defined by equation (25) 

exponent in grid stretch function - see equation (42) 

pressure gradient parameter - see equation (8) 
step size in n coordinate h - An n _-j = n n - h n _| 

5 ' 


boundary- layer thickness measured in the n coordinate. 
Defined as the point where u/U. - .955 unless other- 
wise noted 

boundary- layer thickness measured in the Y coordinate 
nondimensional eddy viscosity - see equation (10a) 
defined by equation (44) 

transformed normal coordinate - see equation (5) 
mixing length - see equation { 1 Od ) 

molecular viscosity 

"k "k 'k 

kinematic viscosity v = p /p 

transformed surface coordinate - see equation (5) 

density 

, , * . * * *, 
wall shear r, = (p 9u /3y ) 

w y = 0 

dummy integration variable - see equation (5) 
coefficients used to compute s-, - see equation (40) 


Subscripts 


e 

w 

00 

n 

exact 


invlscid flow conditions at y = 0 
viscous flow conditions at y = 0 
free-stream conditions 

grid index in n coordinate for finite-difference 
formulation 

value obtained with 640 intervals across boundary layer 


Superscripts 

i present iteration 

1 differentiation with respect to the n coordinate 


- time-averaged quantity 

* dimensional quantity 

-1 inverse of matrix 

-> vector 

T transpose of a matrix (or vector) 

(1),(2), or (3) coefficients associated with unknowns s, f, or v, 

respectively, in equations (32) to (39) 


GOVERNING EQUATIONS 


The method Is demonstrated with application to the steady, two- 
dimensional, incompressible, laminar or time-averaged turbulent boundary 
layer equations^ which are expressed here in Gortler variables^ 2 . 


(* In) " v 2 - 1 > - 2?f || = 0 (momentum) (1 ) 


|~ + 2 ^||’ + T = 0 (continuity) (2) 

with boundary conditions 

f = v = 0 at n 3 0 (no slip, no injection) (3) 


and 


f 1 as n + 00 . 

Other quantities are given by 


(4) 


f* d (x) 

5 = I U G (4») d<j>> n = Y, 


V2? 


(5) 


f = u/U e , 


( 6 ) 


v 

u e 


v + 2 S f l?/ U e> 


(7) 


8 U e dC ■ 


( 8 ) 


Jl 5 1 + E. 


(9) 


The eddy viscosity e is given by 


13 



00 a) 


where 


t * ulfr zi 

( 10 b) 

F - V- exp (-Y + /26) 

( 10 c) 


and 

.085 tanh ( 7^5 -) (k = .41) (lOd) 

The variables in equations (1) to (10) are nondimensional as shown in the list 
of symbols, 

FINITE-DIFFERENCE EQUATION FORMULATION 

The usual approach for reducing a set of differential equations to a set 
of finite-difference equations is to express the partial derivative terms as 
finite-difference quotients and substitute those into the differential equa- 
tions to obtain a set of finite-difference equations. Here we reverse the 
procedure and substitute the differential equations into a finite-difference 
expression to obtain a set of finite-difference equations. I 
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The finite-difference expression around which the method is formulated 
is given by 


9 n - 9 n-l - \ [% * Vl) 



( 11 ) 


where h is the variable step size, n« - n„ i» in the normal coordinate n 
i n m“ I 

and ( ) = 3 ( )/3n with g representing any sufficiently differentiable 

quantity. As mentioned previously, equation (11) is not new. It is the 

fourth-order member of a general class of high-accuracy, two-point methods 

for boundary-value problems discussed by Keller^ 0 . Liniger and Willoughby 14 

analyzed the stability of equation (11) as an implicit method for solving 

initial-value problems for stiff systems of ordinary differential equations. 

2 

Hirsh used the expression to formulate boundary conditions for the three- 
point “compact" scheme applied to the incompressible driven cavity problem. 

A key step in the present procedure, as with the second-order Keller Box 
Scheme, is to reformulate the problem in terms of a first-order set of partial 
differential equations. This is done by defining s = 3f/3n and rewriting 
equations (1) and (2) as 

2 

^ Us - vf) - 2£ ||- + (1 + 6) f 2 - $ (momentum) (12) 



(shear) (13) 


oV _ or it 

~ 35 



If we now define g as a vector 



(continuity) (14) 


( 15 ) 

(16) 
(17) 


10 


then g' will be the right-hand side of equations (12) to (14) and g" its 
derivative with respect to the n coordinate, as follows: 


25 + (1 + 8 ) f 2 


- 3 


9 = 


- 25 - f 

<)5 


(18) 

(19) 

( 20 ) 


2 [25 + (1 + 8 ) (fs) 


g *'(2* 

I 


1) 


-1 


vs + 3(f 2 - 1 ) * 5 - t' |s|s 


- 25 |f-s 


( 21 ) 

( 22 ) 

(23) 


where equation ( 1 ) has been used to express s' in terms of the dependent 
variables f, s, and v (eq. (22)). Substitution of equations (15) to (23) 
into equation ( 11 ) leads to three nonlinear finite-difference equations. 

The 9( )/95 terms in these equations are written so as to handle either 
the Crank-Nicolson scheme, where they are approximated by a central difference 
quotient, or a scheme where they are approximated by three points; thus in 
either case making the equations second-order accurate in the 5 coordinate 
and fourth-order accurate in the normal coordinate n. These nonlinear 
equations are then linearized by Newton's method with the exception of the 
term t'|s|s in equation ( 22 ) and the quantities 6 and s-j which are 
used to compute X and Y + . 

Application of Newton's method yields three linear finite-difference 
equations which can be expressed as 


« _ -b- ->* 

A„z„ 1 + B„z„ = w„ 
n n-l n n n 


(24) 


11 


where 




(25) 


“n = ‘VV/ l 26 > 



a u 

a 12 

a l 3 " 


"hi 

cr 

ro 

b 13 

A n = 

a 21 

a 22 

a 23 

V 

b 21 

b 22 

b 23 


a 31 

a 32 

a 33„ 

1 

_ b 31 

b 32 

b 33_ 


The superscript i in equation (25) denotes the present iteration value. 
A n> B n , and w p are functions of the dependent variables evaluated at 
the i-1 iteration and/or at the previous £- stations. Equations (24) 
must be solved repeatedly until an acceptable level of convergence is 
obtained. Although the double-subscripted matrix elements change values 
with the index n, the index n is suppressed here for simplicity. The 
superscript i will be suppressed from here on for the same reason. 
Appendix A shows how a 21> a 2?) a 23> b^, b 22 , b 23 , and q n are obtained. 

The boundary conditions are 


f-j = v-j - 0 (No slip, no injection) (28) 
and 

f N =1 (29) 

It is noteworthy that, since equations (24) can be applied at the outer 
boundary, the total number of finite-difference equations available exactly 
balances the total number of unknowns. In contrast, a three-point method 
formulated in terms of f, s, and v (since it could not be applied at the 
outer boundary) would result in more unknowns than equations and hence the 


method would require additional boundary conditions or use of three-point 
backward differencing at the outer boundary or a combination of both. 

The resulting linear system of finite-difference equations (24) can be | 

solved in the relatively simple manner which is presented here. 

If equations (24) are applied at n = 2 (n = 1 being the surface and 
n = N locating the outer boundary) and the boundary conditions f-j - v-| = 0 
are applied, three equations with four unknowns are obtained. Thus, we may 
solve for three of the unknowns in terms of the fourth unknown. 

Since f^ is known at the outer boundary, it is convenient to solve for 
s-|, Sg. and v 2 in terms of f 2 , in anticipation of the general form of the 
recursion relations to be derived subsequently. Thus, at n = 2, equa- 
tions (24) can be rewritten as: 


where 


(s 2 , s i , v 2 ) T ~ C 


*2 “ ( b l2 ,b 22 ,b 32^2 


(30 


fb 


C = 


11 


'21 


5 11 


‘21 


'13 


'23 


(31 


'31 


‘31 


'33 


.If we next apply equations (24) at n = 3 and use equations (30) to 
eliminate s 2 and v 2 as unknowns, the result is again three equations and 
four unknowns; namely, f 2 , f 3 , v 3 , and s 3 . Thus, we may solve for three of 
the unknowns In terms of the fourth unknown. Here we solve for f 2 , v 3 , $ 3 , 
in terms of f 3 ; this order being dictated by the choice at n = 2 which was 
made with the consideration of the outer boundary condition f H = 1. If this 
procedure is repeated for n = 4, ..., N, the following general form is 
obtained for the three unknowns written in terms of the fourth: 


. Jk L Jv . 
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M 


= d 0) +e 0> f 


n n 


n n 


(32) 


f , = d J 2) + e < 2) f, 
n-1 n n .1 


n = 3 > 4 » • * • » N • 


(33) 


= d „ (3) + ej 3) f 


n n 


n n 


(34) 


where the superscripts 1, 2, and 3 in parentheses identify the coefficients 
associated with the unknowns s, f, and v, respectively. 

The coefficients d n ^\ d^ 2 ), d„^ 3 ^, e n ^, e„^ 2 ^, and e„^ are 
given by 


n ’ n 


( d n 0> - d n (2)d n (3) ) T “ D iT' [(Pn-'in • r n> T ‘ < a lV a 21> a 31> Td n-) 

- ( a 13" a 23> a 33) Td !i-lJ 


( e n 0) ’ e n (2) ' e n <3> ) T = ' D n'’ (b X2* b 22* b 32) T 


wi th 


V = 


n 


a + a ^ + a e^ 3 ^ 
a 12 + a llVl a 13 e n-l 


M3 


21 


a 22 * a 21 e n-l + a 23 e n-l 


'23 


31 


a 32 + a 31 e n-l + a 33 e n-l b 33 _ 


14 


(35) 


(36) 


(37) 


• 'f 
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i.&fi+'i'i 


Examination of equations (32) to (37) reveals that while the solution for 
the dependent variables s p , f n , and v p cannot be computed as we proceed 

toward the outer boundary, the coefficients d p ^, d p ^, d n ^ 3 ), e n ^, 
e and e ^ can be computed if the initial values dg^), d^ 3 ), 

egO), and e 2 ( 3 ) are known. These initial values are found by comparing 
equations (32) and (34) evaluated at n - 2 with equations (30). The 
initial values are: 



where ip, w are coefficients to compute s-| from the following equation 

s-j = ip + wfg. (40) 

When the outer boundary is reached, equations (32) to (34) can now 
be solved for s^, f ^ , and v^ since the outer boundary condition 

f^ = 1 can be applied. Knowing f M _^, we can next compute s^-j, f N _ 2 , anc * 

v ^-i an(1 likewise all the values of s n , f n , and v proceeding from the 
outer boundary toward the solid boundary. 

RESULTS AND DISCUSSION 

Figure 1 shows the percent error in the wall shear versus the number of 
intervals across the boundary layer, N-l, for the present method, the Spline 4 
method of Rubin and Khosla 6 and the KBS for stagnation point flow (B - 1). The 

number of intervals were 5, 10, 16, 20, 32, 40, 64, 80, 128, 160, and 320. The 

Spline 4 method was unstable for five intervals and thus no point is shown. For 


15 


this case the governing equations reduce to a set of ordinary differential 
equations. The error in the wall shear is defined as 

E x = 100 |s exac t " s ll/ s exact < 41) 

where s exac t is the value of s-j from the higher-order methods with 640 
intervals across the boundary layer. The mesh size for this case varied 
across the layer (e.g., An-j = .05 for N = 11 ) with the last computational 
point being at s 24.2538. The equation used to generate the grid is 
given by 

% - c 0 VO + (42) 

where 

c q = n^| (1 t c-j ) (43) 

= (n- 1 )/ (M-l (44) 

and 

c-j = -.4, a - 8.26 (45) 

The two extra equations needed for the Spline 4 method were taken as 

s N = 0 (46) 

and equation (11) applied at n = 2 to the shear equation. 

Figure 1 shows both the present method and the Spline 4 method of Rubin 
fi 

and Khosla to be fourth-order accurate for a variable mesh. However, for 
the same number of nodal points across the boundary layer, the present method 
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is approximately 10 times as accurate as the Spline 4 method for a variable 
mesh size. To judge the efficiency of the different higher-order methods, 
the computer run times for the different methods were compared for three 
high-accuracy solutions to this problem; namely, E x = .15, .015, and .0015. 

To achieve an error of E « .15 'the present method, the Spline 4 method, 
and the KBS solutions were computed with 7, 11, and 27 nodal points, respec- 
tively. (See fig. 1.) The convergence criteria for these three cases required 

the value of E between consecutive iterations to be accurate to two signif- 
x 

icant figures. Results are summarized in tables I, II, and III. These tables 
show the present method to be significantly faster than the Spline 4 and the 
KBS to achieve the same degree of accuracy. It should be noted that a key 
reason for the greater efficiency of the present method over the Spline 4 
method is due to the present equation formulation which allows the resulting 
finite-difference equations to be solved efficiently in a coupled manner using 
Newton's method. 

Figure 2 shows a plot of the percentage error in the wall shear versus 
the number of intervals across the boundary layer for a model turbulent 
problem for the present method and the KBS. The equations for this case were 
obtained from the nonsimilar equations by setting 6=0 and 9( )/35 =0. 

The Reynolds number was 1.88 million, £ = 1, 6* = 24.5 with the outer boundary 
located at s 60. Convergent solutions using the Spline 4 could not be 
obtained for this case. (The possibility of a coding error exists in the 
program and this will be investigated further.) The constants in equation (42) 
for the ti grid distribution are c-j s .05 and a - -109. Figure 2 shows the 
present method to be slightly better than third-order accurate for the model 
turbulent problem. In order to determine the efficiency of the present method, 
solutions were computed with both methods for a .1-percent error in the wall 
shear. The present method requires approximately 21 nodal points and the 
KBS 92 nodal points. Results show the present method to be approximately 
3.83 times faster than the KBS to achieve the same order of accuracy. All 
calculations were made on a CDC Cyber 175 computer with a FORTRAN Extended 
Compiler, optimization level = 0. The same calculations would require about 
2.5 times more CPU time on a CDC 6600. No under-relaxation was required for 
either test case. 
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CONCLUDING REMARKS 


A fourth-order box method for calculating high-accuracy numerical solu- 
tions to parabolic, partial differential equations in two variables or ordinary 
differential equations has been presented. The method is the natural extension 
of the second-order Keller "Box" scheme to fourth order and has been demon- 
strated with application to the incompressible, laminar and turbulent boundary- 
layer equations. Numerical results for high accuracy test cases show the 
present method to be significantly faster than other higher-order and second- 
order methods available today to achieve the same accuracy. For the test 
cases reported on the governing equations reduced to ordinary differential 
equations. Calculations for nonsimilar cases have been completed and will 
be included in a later report. 

Richardson extrapolation to zero mesh size has not been considered in 
the present paper. Keller has shown this technique to be a valuable one 
for improving the accuracy of the KBS; hence, a final comparison must be 
made on this basis. 
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APPENDIX A 




SHEAR EQUATION COEFFICIENTS 


We show here how to identify the matrix elements of A n > B n and the 
vector component of w n for equations (24). The procedure is illustrated for 
the shear equation for laminar, similar flow (A = 1, 3( )/3£; * 0) where for 
the shear equation 


g - f 


g = s 


g" = s' « vs + B (f 2 - 1) 


Substitution of equations (Al) to (A3) into equation (11) evaluated at the 
ith iteration yields the following nonlinear finite-difference equation 


- K.i -e|( f2 Ci -ij] = 


Equation (A4) is now linearized via Newton's method which can be shown to 

4 p \ 

correspond to the following linearization for the products (vs) and (f ) ; 


(vs )^ = v^s 1 + vV " 1 - (vs ) 1 - 1 


21 


.W/, 


c? a •• Jiy ; *.vs W- ■ W-I» s.v r -■*■: -4 




(f 2 ) 1 = zf^V - (f 2 ) 1 " 1 


*A6) 


Applying equations (A5) and (A6) to equation (A4), the linear finite-difference 
equation can be written as 

a 21 s n-1 + a 22 f n-l + a 23 v n-l + b 21 s n + b 22 f n + b 23 v n = q n ^ A7) 


where 


‘21 


K/9 h 2 v,i‘l 

h/2 • 12 Vi 


(AS) 


‘22 


1 . lhlfi-1 
1 6 f n-l 


(A9) 


1 - hi c i"l 

*23 12 s n-l 


(A10) 


J 21 


h + si J-i 

2 12 v n 


(All) 


J 22 


1 + fi"l 
1 6 T n 


(A12) 


23 


J-l 
12 s n 


(A13) 


and 


„ -i. 

q n ' 12 




(A14) 
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Figure 2.- Percent error in the wall shear versus the number 
of intervals across the boundary layer. 




